The merger of white dwarf-neutron star binaries: 
Prelude to hydrodynamic simulations in general relativity 
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White dwarf-neutron star binaries generate detectable gravitational radiation. We construct 
Newtonian equilibrium models of corotational white dwarf-neutron star (WDNS) binaries in circu- 
lar orbit and find that these models terminate at the Roche limit. At this point the binary will 
undergo either stable mass transfer (SMT) and evolve on a secular timescale, or unstable mass 
transfer (UMT), which results in the tidal disruption of the WD. The path a given binary will 
follow depends primarily on its mass ratio. We analyze the fate of known WDNS binaries and 
use population synthesis results to estimate the number of LISA-resolved galactic binaries that will 
undergo either SMT or UMT. We model the quasistationary SMT epoch by solving a set of simple 
ordinary differential equations and compute the corresponding gravitational waveforms. Finally, we 
discuss in general terms the possible fate of binaries that undergo UMT and construct approximate 
Newtonian equilibrium configurations of merged WDNS remnants. We use these configurations to 
assess plausible outcomes of our future, fully relativistic simulations of these systems. If sufficient 
WD debris lands on the NS, the remnant may collapse, whereby the gravitational waves from the 
inspiral, merger, and collapse phases will sweep from LISA through LIGO frequency bands. If the 
debris forms a disk about the NS, it may fragment and form planets. 

PACS numbers: 04.25.D-,04.25.dk,04.30.-w 



I. INTRODUCTION 

The inspiral and merger of compact binaries represent 
some of the most promising sources of gravitational waves 
(GWs) for detection by ground-based laser interferome- 
ters like LIGO 0,0, VIRGO GEO H, TAMA 0,0 
and AIGO [|[, as well as by proposed space-based inter- 
ferometers like LISA [|] and DECIGO [Io|. Extracting 
physical information from gravitational radiation emitted 
by compact binaries requires careful modeling of these 
systems (see [ll| for a review). Most effort to date has 
focused on modeling black hole-black hole_(BHBH) bina- 
ries [H [H Q [S Hi E3, El d, H ID, MM_ and neu- 
tron star-neutron star (NSNS) binaries [Mill, [H, H3, 
HI, IH H3] , with some recent work on black hole-neutron 
star (BHNS) binaries 0, H H H M, M, [13, M, HI . 

In this paper we begin to explore WDNS binaries. 
They are plausible sources of low-frequency GWs for 
LISA, DECIGO and possibly also, as we shall see, high- 
frequency GWs for LIGO, VIRGO, GEO, TAMA and 
AIGO. 

Like NSNS binaries, WDNS binaries are known to ex- 
ist. We have compiled a list of observed WDNS bina- 
ries in Tables Q] and [Til Table [TJ lists those binaries with 
relatively well-determined individual component masses. 
By "relatively well-determined" we mean that in deter- 
mining the masses no assumption was made about the 
mass of the NS. Table ILTI lists those binaries for which the 
masses of the individual components have not been well- 
determined as yet. For these cases the NSs have been as- 
sumed to have a mass of either 1.35Af or 1.40M Q . Of all 
the objects presented in Tables HI and UTI only B2303+46, 



J1141-6545, J0751+1807 and J1757-5322 have been pos- 
itively identified as a pulsar with a WD companion. For 
all other systems there is at best strong evidence that the 
companion of the observed pulsar is a WD. 

The frequency of GWs emitted by the binary with the 
shortest period, J1T41-6545, is ~ 1.2 x 10~ 4 Hz. This 
frequency lies in the expected band of LISA, but well be- 
low the cutoff frequency of ~ 10 _3 Hz due to the double 
WD confusion background [I(| • It would therefore be im- 
possible to resolve this signal. Assuming the quadrupole 
approximation for a WDNS binary with masses Mwdj 
Mns, reduced mass fi, and total mass Mr = A7ns + Mwd 
in circular orbit with angular frequency f2, the amplitude 
of GWs coming from this object is h = 4/xMt / '(rA) [4l| . 
where r is the distance to the object and A the binary sep- 
aration. In [42] a lower bound is given to the distance to 
PSR J1141-6545 of r > 3.7 kpc. Using the data of Table 
UJ we find that the maximum amplitude of the expected 
waves coming from this object is h max = 1.36 x 1CU 23 . A 
gravitational wave of this amplitude is below LISA's sen- 
sitivity at 10 _4 Hz (h s ~ 10~ 21 ) and hence undetectable. 
In order for detectable gravitational wave signals to be 
emitted the orbital separation of PSR J1141— 6545 has 
to decrease by a factor of about 100. 

The emission of gravitational radiation will cause a 
binary to inspiral to a close, nearly circular orbit. In 
this regime the binary can be described by a sequence of 
quasiequilibrium configurations. Once the separation of 
a WDNS binary becomes a few times larger than the size 
of the WD, the amplitude of the emitted GWs, as we will 
see, is large enough to allow for detection by space-based 
gravitational wave observatories. An important issue to 
address then is the number of Galactic WDNS binaries 
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that could be detected by space-based interferometers. 

Population synthesis calculations by Nelemans et al. 
[H show that there are about 2.2 x 10 6 WDNS bina- 
ries in our Galaxy, and that they have a merger rate of 
1.4 x 10 _4 yr -1 . Furthermore, Nelemans et al. find that 
after a year of integration, LISA should be able to detect 
128 WDNS binaries and, after considering the contribu- 
tion of the double WD background GW noise, resolve 38 
of these. On the other hand, calculations by Cooray [50| . 
give much more conservative numbers of resolved WDNS 
binaries. In particular, Cooray finds that the number of 
LISA-resolved WDNS binaries ranges between 1-10, us- 
ing a WDNS merger rate between 10 _6 yr _1 -10 _5 yr _1 . 
Cooray's upper limit was based on merger rates calcu- 
lated by Kim et al. [Fll ]. 

In this work we focus on WDNS binaries in close bi- 
nary separations, and examine the termination point of 
quasiequilibrium sequences describing such binaries. Sev- 
eral different astrophysical scenarios can result in such a 
termination point for binaries in general, including direct 
contact of the binary components, Roche lobe overflow 
by one of the two companions, or the binary reaching an 
innermost stable circular orbit (ISCO). As we will find, 
quasiequilibrium sequences of WDNS terminate when the 
WD fills its Roche lobe, resulting in mass transfer from 
the WD onto the NS across the inner Lagrange point. 
This mass transfer can either be stable (SMT) or unsta- 
ble (UMT). We also refer to the latter scenario as the 
tidal disruption of the WD by the NS. 

To determine which of these two outcomes is likely for a 
given system, we will follow the approach of Verbunt and 
Rappaport [52| and Faber et al. 53 1. As indicated in our 
Tables Q] and HH we shall find that among the observed 
WDNS binaries, some will undergo SMT, while others 
will undergo tidal disruption (i.e. UMT). Note that mass 
transfer stability has also been studied in [54l . |55| . 

In the case of SMT, the orbital evolution of the binary 
occurs on a secular timescale determined by the emis- 
sion of gravitational radiation. Therefore, the quasista- 
tionary conservative treatment of Clark and Eardley (5(| 
and Faber et al. [53| is adequate to follow the evolution 
during this secular phase. Note that a quasistationary 
treatment to follow the orbital evolution of binary sys- 
tems has been employed by several authors in the past. 
For example, Rappaport et al. [13] studied compact bina- 
ries where the mass of the secondary can be up to ~ 1-Mg 
and modeled as an — 3/2 polytrope. Fryer et al. [58| 
employed a non-conservative quasistationary approach to 
study the evolution of white dwarf-black hole binaries, 
while Marsh et al. [55J studied the evolution of double 
WD binary systems. Both the tidal disruption and SMT 
)hase of double WD systems have also been studied in 
M H M, HI H HI via SPH simulations and in 
6g, |67[ via grid-based hydrodynamic calculations, all in 
Newtonian gravitation. Finally, Newtonian SPH simula- 
tions of encounters of WDs with intermediate mass BHs 
(Mbh ~ 10 3 M©) have been performed in (68l. [69|. 

In the case of tidal disruption, on the other hand, 



the system will evolve on a hydrodynamical (orbital) 
timescale. In this scenario the NS may plunge into the 
WD and spiral toward the center of the star liberat- 
ing its gravitational potential energy as heat in the WD 
material. Alternatively, the NS may be the receptacle 
of massive debris from the disrupted WD. Depending 
on the details of the equation of state, a cold degener- 
ate gas can support a maximum NS rest-mass between 
1 .89AT Q — 2.67M (corresponding to a gravitational mass 
between 1.65M© — 2.20M Q ) against catastrophic collapse 
if it is not rotating (the OV limit), about 20% more mass 
if it is rotating uniformly (a "supramassive NS" , e.g. [z3]), 
and at least 50% more mass if it rotates differentially (a 
"hypermassive NS") [zE ff^. The fate of the merged 
WDNS then depends on the initial masses of the pro- 
genitor stars, the degree of mass and angular momentum 
loss during the WD disruption and binary merger phases, 
the angular momentum profile of the NS remnant and the 
extent to which the remnant gas is heated by shocks as 
it pours onto the NS and forms an extended, massive 
mantle. These are issues that require a hydrodynamic 
simulation to resolve. 

Moreover, ascertaining whether or not the neutron star 
ultimately undergoes a catastrophic collapse to a black 
hole (either prompt or delayed) requires that such a sim- 
ulation be performed in full general relativity. We plan 
to explore some of these alternative hydrodynamical sce- 
narios in detail in the future, aided by simulations that 
employ our adaptive mesh refinement (AMR) relativistic 
hydrodynamics code [38[. 

In this paper we survey the problem in qualitative 
terms. In Section [TT] we identify some of physical param- 
eters that are likely to play a key role in the evolution 
of a WDNS binary system. We then model the secu- 
lar inspiral epoch of the binary by constructing Newto- 
nian equilibrium models of corotational WDNS binaries 
in close circular orbit, up to the Roche limit. In Sec- 
tion IIIII we follow the stability analysis of Verbunt and 
Rappaport 52] in order to determine the late-evolution 
of these binaries (SMT vs. tidal disruption). We treat the 
quasistationary SM T ep och by applying the approach of 
Clark and Eardley [56j and Faber et al. [53| and com- 
pute the corresponding gravitational waveforms. We dis- 
cuss in general terms the possible outcomes of binaries 
that undergo UMT and construct approximate equilib- 
rium configurations of merged WDNS remnants. We use 
these models to make some predictions regarding our fu- 
ture fully relativistic simulations of these systems. Fi- 
nally, we conclude in Section IIV1 where we summarize 
the main findings of this work. 



II. EQUILIBRIUM CONFIGURATIONS 

In this section we follow Faber et al. [Hj], who stud- 
ied BHNS systems, to distinguish the different timescales 
which are relevant for WDNS systems. In particular, we 
will find that WDs are likely to rotate synchronously with 
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TABLE I: WDNS binaries with well-determined masses. From left to right the columns give the name of the object, the orbital 
period, the corresponding quadrupole GW frequency, the WD mass, the NS mass, the total mass, the mass ratio q — A/wd/A^ns 
and the stability of mass transfer after its onset. 



PSR 


T (days) 


/gw = f (10" 4 Hz) 


Mwd(M q ) 


M NS (M ) 


M T (M©) 


q 


Stable? 


B2303+46 a 


12.34 


0.0187 


1.3 


1.34 


2.64 


0.97 


No 


J0621+1002 6 


8.32 


0.0278 


0.67 


1.70 


2.37 


0.394 


? 


J1141-6545 c ' d 


0.198 


1.169 


1.02 


1.27 


2.29 


0.803 


No 


B1516+02B 6 


6.858 


0.0337 


0.13 


2.08 


2.21 


0.0625 


Yes 


J1713+0747" 


67.8 


0.0034 


0.33 


1.60 


1.93 


0.206 


? 


B1855+09 a 


12.3 


0.0188 


0.267 


1.58 


1.847 


0.169 


Yes 


J0437-4715 11 


5.74 


0.0403 


0.236 


1.58 


1.816 


0.149 


Yes 


■11012+5307° 


0.605 


0.382 


0.16 


1.64 


1.80 


0.097 


Yes 


■10751+1807° 


0.263 


0.88 


0.125 


1.26 


1.385 


0.099 


Yes 



"Stairs |4."J and references therein. 
6 Nice et al. H. 
c Bhat et al. 1431 
d Bailes et al. [491 . 
Treire et al. [4711 . 



TABLE II: WDNS binaries which do not have well-determined masses. From left to right the columns give the name of the 
object, the orbital period, the corresponding quadrupole GW frequency, the WD mass, the NS mass, the total mass, the mass 
ratio q = Mwd/A/ns and the stability of mass transfer after its onset. 



PSR 


T (days) 


/gw = f (10- 4 Hz) 


Mwd(M q ) 


Mns(M ) 


Mt(Mq) 


<7 


Stable? 


.11435-60° 


1.355 


0.1708 


1.10 


1.40 


2.50 


0.785 


No 


■11157-5114* 


3.507 


0.066 


1.14 


1.35 


2.49 


0.844 


No 


.11453-58° 


12.42 


0.0186 


1.07 


1.40 


2.47 


0.764 


No 


J1022+1001° 


7.805 


0.0296 


0.872 


1.40 


2.272 


0.623 


No 


B0655+64 a 


1.029 


0.2948 


0.814 


1.40 


2.214 


0.581 


No 


J2145-0750° 


6.839 


0.0338 


0.515 


1.40 


1.915 


0.368 


7 


J1757-5322 6 


0.453 


0.511 


0.55 


1.35 


1.90 


0.407 


7 


J1603-7202° 


6.309 


0.0367 


0.346 


1.40 


1.746 


0.247 


7 


J1810-2005° 


15.01 


0.0154 


0.34 


1.40 


1.74 


0.243 


7 


J1904+04° 


15.75 


0.0147 


0.27 


1.40 


1.67 


0.193 


7 


J1232-6501° 


3.507 


0.066 


0.175 


1.40 


1.575 


0.125 


Yes 



"Tauris et al. |4S| and references therein. 
'Edwards and Bailes [4gj . 

compact companions. We then construct and describe 
the basic features of Newtonian equilibrium configura- 
tions of corotational WDNS binaries in circular orbit. 



A. Timescales 

We begin by defining the primary as the NS and the 
secondary as the WD. The relevant timescales are the 
orbital period T, the gravitational wave timescale of the 
binary tow an d the viscous timescale i v ; s of the WD. 
We use the dynamical timescale idyn of the secondary in 



order to rescale all other timescales. Before we proceed 
we define the mass ratio q of the binary according to 

1 = — 7 — • (1) 

Mns 

We also define the compaction parameter C of the WD, 

C=^R, (2) 

where i?wD is the WD radius. 

By virtue of equations (fTJ) and @ we can define the 
dynamical timescale of the WD (i.e. the time it takes a 
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0.2 0.4 0.6 0.8 1 1.2 1.4 0.2 0.4 0.6 0.8 1 

M WD /M q 

(a) (b) 



FIG. 1: (a) Compaction C = Mwd/Rwd of isolated WD equilibrium models plotted as a function of the WD mass. The 
adopted EOS is that of a cold degenerate ideal electron gas for fi e — 2. (b) Critical turbulent viscosity parameter for a WD 
plotted against the mass ratio of a WDNS binary for various values of the WD compaction. 



sound wave to propagate across the star « the free-fall 
timescale) as 



tdy 



^WD 



M 



C- 3/2 M WD = qC- 3/2 M } 



NS- 



WD 



(3) 



Note that we use geometrized units (G = c = 1) here 
and throughout this paper with few exceptions where we 
explicitly restore cgs units. 

We are interested in WDNS systems that have reached 
their Roche limit. Therefore, we need an expression for 
the Roche lobe "radius" R Y of the WD. There are two 
approximations to the Roche lobe that have been used 
widely in the WD literature, namely the Paczyhski [73| 
and Eggleton [74[ formulae. We can write these in the 
general form 



Rr=Af(q), 



(4) 



where A is the orbital separation. In the Paczyhski ap- 
proximation we have 



/(<?) = co 



Q 



1 + 9 



1/3 



(5) 



with Co = 0.46, whereas in the Eggleton approximation 
we have 



/(<?) = 



ciq 



2/3 



c 2 g 2 / 3 + ln(l + g 1 /3)' 



(6) 



with a = 0.49 and c 2 = 0.6. 

Clearly, the Paczyhski approximation is simpler than 
the Eggleton approximation, but the latter is more accu- 
rate than the former [74J . For the purposes of this section 
it is adequate to use the Paczyhski approximation (0, 



but for our treatment of mass transfer in the following 
section we will adopt the Eggleton approximation ([6]). 

The WD reaches its Roche limit when the Roche lobe 
radius becomes equal to the WD radius, R T = i?wD- This 
occurs at orbital separation An given by 



A R = csV 1/3 (l + q^C^Mwd, 
= Co - 1 g 2 / 3 (l + (? ) 1 / 3 C- 1 M NS . 

At this point the Keplerian orbital period is given by 



(7) 



T= ^ = 2n yik = 2 ^o V2 qC-^M NS , (8) 

where is the orbital angular frequency and Mp = 
Mwd + Mns the total mass of the system, so that 



T 

^dyn 



2-kc, 



-3/2 



20.1 



(9) 



Treating the components of the binary as point masses 
and assuming the quadrupole approximation, we obtain 
the gravitational wave timescale 



.4 



A 4 



(10) 



A 64 MnsMwdMt ' 

where an overdot stands for a time derivative. At the crit- 
ical (Roche limit) orbital separation the GW timescale is 

5 „ 

(11) 



64 M ns M W dM t : 



^c -V /3 (i + g) 1/3 cr 4 M NS , 



and hence 



^ = A c -V /3 d + 9 ) 1/3 c- 5 / 2 

idyn t>4 



(12) 
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For a realistic (massive) WD of about a solar mass, radius 
of 10 3 km (i.e. C ~ 10~ 3 ) and a NS companion such 
that q £ [0.3, 1] the above ratio is of order 10 7 . Thus, 
the GW timescale is about 10 6 times the orbital period 
(see Eq. (|9])). Therefore, during the inspiral epoch (up 
until reaching the Roche limit) the evolution occurs on a 
secular timescale and the binary system can be taken to 
be in a quasiequilibrium state. 

We now turn to the viscous timescale i v ; s of the 
WD. Viscosity is particularly important in determining 
whether the WD will be rotating synchronously in the 
WDNS binary or not. According to [H, [75[, if the vis- 
cous timescale is short enough, tidal dissipation can syn- 
chronize the binary. In particular, this will occur if 



Rwd 



> 60(l + 9 ) 5/ V 2/3 C* 3 , 



(13) 



or 



l(l + (7 )- 5 /V/ 3 C- 4 M NS . (14) 
60 



According to Horedt [76[, turbulence is very likely to be 
tidally induced on the WD by the NS, when the difference 
between the orbital angular velocity £1 and the spin an- 
gular velocity of the WD u> is such that \Q — lj\ > 0.157. If 
we then assume a-law turbulent viscosity as the primary 
source of viscosity, so that 



«vis 



^dyn 
f . ' 

< VIS 



(15) 



then synchronization occurs when 

a vis > (Writ = 60(1 + q) 5 ^q- 2 / 3 C^ 2 . (16) 

Therefore, for a given WD compaction C, a v is,crit in- 
creases with decreasing q. 

The compaction of stable WDs supported by cold 
degenerate ideal electrons and obeying the Tolman- 
Oppenheimer-Volkov (TOV) equations [4l[ lies in the in- 
terval [10 -6 , 10~ 3 ] (see Fig. 1(a) here and in the follow- 



ing we assume a mean molecular weight per electron of 
He = 2). In Figure 1(b) we plot the critical turbulent 
viscosity parameter a V is,crit versus q for various values of 
the WD compaction. For WD compactions C < 5 • 10~ 4 , 
which implies A/wd 1.3A/©, very small turbulent vis- 
cosity (a v is — 10~ 6 ) is adequate to satisfy condition 
(fTB)) and hence tidally lock the WD. Plausible values 
of turbulent viscosity lie in the range a v ; s £ [0.001,0.1] 
fn\ . Therefore, we conclude that the WD component 
of WDNS binaries will most likely be corotational. It is 
interesting to contrast this result to the case of a NS in 
a BHNS binary, for which Faber et al. [53| found more 
likely to be irrotational. 

Purely hydrodynamic turbulence may not be the only 
source of viscosity for WDs. Magnetic fields may drive 
MHD turbulence and thereby lock the WD tidally as well. 



B. Basic equilibrium equations 

We now describe the construction of corotational equi- 
librium models of WDNS binary systems in order to 
model the early quasistationary inspiral epoch of the bi- 
nary. 



1. Assumptions 

Here we adopt Newtonian gravitation to model the 
equilibrium tidal distortion of the WD. This greatly sim- 
plifies the formulation, and is reasonable due to the low 
compaction of WDs, C < 10" 3 . The leading order rela- 
tivistic corrections will be of this order, and we therefore 
choose to neglect them. 

To further simplify the analysis we model the NS as a 
point mass. The size of the NS is negligible compared to 
the much larger WD, i?NS -Rwd, and plays little role in 
determining the WD profile. Given the large compaction 
of the NS (Cns ~ 0.1), its deviation from sphericity will 
always remain very small. 





10 _o 




10" 10 




10" 12 


E 


10" 14 


_^ 


10" 16 


Q- 


10" 18 




10" 20 




10" 22 




10" 24 



EOS 
r=4/3 polytrope 
r=5/3 polytrope 



10 -16 10 -14 10 -12 10 -10 10 _1 
p/u„ (1/km 2 ) 



FIG. 2: The exact ideal degenerate electron EOS approaches 
a r = 5/3 polytrope in its low density (nonrelativistic elec- 
tron) limit, and a F = 4/3 polytrope in the high den- 
sity (ultra-relativistic) regime. The transition occurs at 
p/n? « l(T 12 km" 2 = 1.347 x 10 6 g ■ cm" 3 . 



2. Equation of State 

We model the WD as a synchronously rotating self- 
gravitating fluid supported by degenerate electron gas 
pressure. We assume zero temperature and ignore elec- 
trostatic corrections and effects of inverse /3-decay and 
model the degenerate electrons as an ideal Fermi gas. We 
assume a mean molecular weight per electron of \i e = 2, 
as would be the case for a typical C-0 star. While the 
pressure is generated entirely by the degenerate electrons, 
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the rest mass is contributed by the baryons. In the ex- 
treme low and high density limits, corresponding to the 
nonrelativistic and the ultrarelativistic electron regime 
respectively, this EOS state can be approximated by a 
polytropic relation 

P = Kp r , (17) 

where P is the pressure, p is the (rest-mass) density, and 
K and T are constants. We denote the corresponding 
internal energy density ej. We graph the exact EOS, 
together with its two asymptotic limits, in Figure[2l Note 
that restoring cgs units lKm~ 2 = 1.347 x 10 18 g • cm 3 . 



1. Implementation of EOS 

We implemented the EOS in tabular form, using poly- 
nomial interpolation between tabulated values. The en- 
thalpy is computed using the tabulated values of p and 
P by integrating Eq. We adopted a fourth order 

Runge-Kutta method [8l| to carry out this integration. 
The enthalpy h is then tabulated in the final EOS table 
which our code reads at run time. 



2. Coordinate set up 



3. Equilibrium Equations 

An equilibrium configuration of a WDNS binary is de- 
termined by the Poisson equation, 



V 2 = 47T/9, 



(18) 



and the integrated Euler equation, which, assuming coro- 
tation, can be written 



h + <l>- 1 -tfrl t =C. 



(19) 



Here r rot is the distance from the binary's axis of rotation 
to a given fluid element, C is a constant, and h is the 
specific enthalpy defined as 



J P 

We can also write Eq. (fT9"|) as 

h=C -<f> eS 

where 

1 

2 



S = (f>- ^ 2 r r 2 ot 



(20) 

(21) 
(22) 



is the effective potential. 

The gravitational potential of the point mass NS is 



0NS 



Mi 



NS 



(23) 



where tns is the distance from the NS. The gravitational 
potential (f> at any given point then can be written 



t>WD + 0NS = 0WD 



Mm 



(24) 



where </>wd is the potential due to the WD. 



C. Numerical Technique 

We construct numerical solutions to Eqs. (fTS)) and (|2"Tj) 
using an algorithm based on [78l Pzoj ] . Our implementa- 
tion is very similar to that described in [80| . except that 
here we adopt a tabulated EOS and allow for binary com- 
panions of similar mass. 



We define Cartesian coordinates in the comoving frame 
of the binary following We assume the binary's axis 
of rotation to be parallel to the z-axis, and located at 
x = x rot and y = (see Figure [3]) . 

It is numerically convenient to have the surface of the 
WD intersect the x-axis at fixed coordinate locations 
xa and xb, where xa is the coordinate of the point 
closer to the axis of rotation. To achieve this we non- 
dimensionalize the coordinates using the x dimension of 
the WD, r e . For a spherical WD, r e is its radius, oth- 
erwise we calculate it as half the WD's diameter in the 
x-direction. We then define the grid coordinates as 



x = x/r e , y = y/r e , z = z/r 



(25) 



We denote grid coordinates with hats, and their physi- 
cal counterparts without hats. We choose xa = — 1 and 
xb = 1 to define the limits of the star in the grid coordi- 
nates. 

We also define 



£>WD 



<hvr>/rl, A = fi/r e 



(26) 



The Laplace operator on the grid, V 2 , scales as V 2 = 
r 2 V 2 . 

Given these definitions and utilizing the fact that the 
potential of the NS is known, we rewrite the Poisson Eq. 
(fl8)l with respect to the computational grid as 



'WD 



(27) 



where p here is the rest-mass density of the WD. We also 
write the effective potential, Eq. (f22|) . as 



t>cs = 4> - 



(x-x ro t) +y 



(28) 



where 



r e (pwD 



Finally, we combine Eqs. 
integrated Euler Eq. (|2"Tj) as 



TU- 



NS 2 1 



-n 2 



(29) 

r e rNS 

and (j2"9"| to rewrite the 



(x-x rot f + y 2 +C. (30) 
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FIG. 3: The coordinate system used in our calculations. We specify that the WD spans from x = — 1 to x = 1 along the 
a;-axis, and rescale with each iteration. The NS is modeled as a point mass on the x-axis at xns- Our computational grid takes 
advantage of symmetries along the y and z axis, and only calculates the positive values here. 



Given that the gravitational potential of the neutron star 
is known analytically, we may restrict the numerical grid 
to a region around the WD (see also (SOj). In our appli- 
cations we impose al/r fall-off condition on 0wd at the 
outer boundaries, located at coordinate values of x = ±4 
and y, z = 4. Note that for large binary separations the 
neutron star is outside the numerical grid. 



3. Iterative solution 

We compute equilibrium models of WDNS binaries by 
simultaneously solving Eqs. (j2"?| and (I30[) using an al- 
gorithm that is described in detail in [S0|] . For a given 
binary separation, mass ratio q and central density of the 
WD p c this iteration also determines the constants C, r e 
and the orbital angular frequency fl. 



4- Constant mass sequences 

In an effort to mimic the inspiral of the binary as it 
looses angular momentum to gravitational radiation, we 
construct sequences of equilibrium models with constant 
WD and NS mass and decreasing orbital separation A. 
This approach is reasonable because for large orbital sep- 
arations the inspiral occurs on a secular timescale, so that 
at any given time, the binary can be taken to be an equi- 
librium configuration. We calculate the separation from 
the point source NS to the center of mass of the WD as, 



A = r e (x N s - iwD.cm)- 



(31) 



The position of the NS on the coordinate grid Xns is 
stepped in sequence from large spacing (where the WD 
is nearly spherical) to the termination of the equilibrium 
sequence. The binary separation is calculated at each 
step from Eq. (f5T|) . 



D. Diagnostics and tests 

Below we summarize several physical diagnostics and 
tests which we use to verify the reliability of our code. 

1. Diagnostics 
We calculate the total energy according to 



E cq = T + W + U. 
In Eq. (f3"2"| T is the kinetic energy, 

T = — I v z dm 



(32) 



- ^ 2 ^ot,NS^NS + \& J r* ot pd 3 x, (33) 

WD 

where r ro t,NS and r ro t are the distances of the NS and a 
WD fluid element from the axis of rotation respectively. 
The gravitational potential energy W is given by, 

=~ J M NS 5(x-x NS )<f>d 3 x + ^ J p<j)d 3 x 



NS 



WD 



iM NS 0wD(x N s) + \ I /°(^NS + (j)WD)d 3 X, 



WD 



(34) 

where </>wd(xns) ~ -Mwd/ins is the WD potential at 
xnS: the position of the NS. The internal energy term U 
in Eq. gl is 



U = J Eid 3 x, 

WD 



(35) 
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FIG. 4: Normalized virial V norm for J1141-6545 showing virial 
equilibrium to well within 1% for the coarsest grid resolution, 
and second order convergence. Three grid resolutions are plot- 
ted, with 64x32x32, 96x48x48, and 128x64x64 grid points. 



where ti is the internal energy density of the fluid. At in- 
finite binary separation, which corresponds to the special 
case of a WD in isolation, Eq. (|3"2")l reduces to 



£■00 = Woo + Uac, 



(36) 



where Woo is the potential energy of a spherical WD. We 
define the system's binding energy Eb as 



Eb — E cq — E a 



(37) 



Similarly, we calculate the angular momentum of the sys- 
tem, J, as 

J = JNS + JWD = ^^NS^rot.NS + ^ / r^pdPx. ( 38 ) 



WD 



2. A convergence test: the virial theorem 

An equilibrium binary should satisfy the virial theo- 
rem. We calculate the virial expression as an independent 
diagnostic test of the accuracy of our code. Following 
[4lj . the virial equation for our binary is 



2T + W + 311 = 0, 

where T and W are given by Eqs. 
pressure term II is 



and 



n 



Pd 3 x. 



(39) 
The 

(40) 



WD 



A normalized virial relation that quantifies the numer- 
ical error in our calculations, is then 



V m 



2T + W + 3IT 



\2T\ + \W\ 



3H 



(41) 



The above virial expression should converge to zero with 
increasing resolution and reflect the second-order accu- 
rate finite difference representation of the Poisson equa- 
tion which our code employs. In Figure [5] we show a con- 
vergence test based on the normalized virial expression 
for one of the constant mass sequences described in Sec- 
tion QTCMl demonstrating that our code is second-order 
convergent. 



3. WDs in isolation 

As another test of our code, we studied WDs in iso- 
lation and compared the results against those computed 
with a ID code that takes advantage of the spherical sym- 
metry of the WD in isolation. To study an isolated WD 
with our 3D code we simply set the mass of the NS to 
zero. 

In Fig. [5] we show the mass and radius of the WD as 
a function of central density. We plot results from our 
3D code, and show that they are in agreement with those 
obtained with our ID code. We label the latter, which 
can be obtained to machine accuracy, as "exact" . In 
Fig. [5] we also show that the masses calculated with our 
3D code converge to second order toward the exact ID 
results. 



4- Orbital Angular Velocity 

In Fig. [7] we graph the orbital angular velocity Q versus 
the binary separation A, for two constant mass sequences 
representing the binaries J1141-6545 and B1855+09 
(cf. Table[T|. The binary J1141-6545 is a representative of 
the high-mass WD population, with relatively equal WD 
and NS masses, while B1855+09 represents the low-mass 
WD population (compare [43l]). 

Figure[7]demonstrates that the orbital angular velocity 
follows Kepler's third law 



n 2 = 



(42) 



very closely, even when the WD is significantly distorted 
at small separations. This is similar to the result of [6lj ] . 
who found that only at very small separations did the 
angular speed of identical binary polytropes (similar to 
a low mass double WD system) measurably diverge from 
Kepler's law. This behavior is a consequence of the WD 
being centrally condensed, so that most of its mass is 
concentrated close to its center. Therefore, the orbital 
motion is well approximated by that for point masses, 
even when the envelope is significantly distorted. 



E. Numerical Results 

With our code we produce equilibrium models of 
the WDNS system and quasiequilibrium sequences of 
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FIG. 5: Mass M (a) and radius R (b) of a spherical WD versus their central density p c , for the EOS described in Section 
HTBI Results from our 3D code at three resolutions (corresponding to 8, 16, and 32 points across the WD radius) are plotted 
together with the "exact" results from the ID code. 
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Roche limit we calculate numerically the volume V of the 
WD and define the volume radius R as 



Pc ( km "' 



FIG. 6: Relative difference between the ID (exact) calculation 
of the mass for spherical WDs and the corrresponding 3D 
results vs. density. Here AM = |M — Msd|, where M is the 
ID result. The plot demonstrates second order convergence 
of the 3D solutions towards those of the exact ID integration. 



constant mass that mimic the evolutionary inspiral se- 
quences. 



1. Roche limit 

At a certain separation our code ceases to converge 
and no equilibrium models can be constructed for smaller 
separations. At this critical separation the WD forms 
a cusp. This critical separation is the familiar Roche 
limit, Ar. It is convenient to normalize the Roche limit 
separation with the WD volume radius, i.e., the radius 
of a sphere which has the same volume as the WD. At 



R 



W 

ill 



1/3 



(43) 



In Table HTT1 we tabulate the critical separation in units 
of the volume radius of the WD as found from the nu- 
merical model and as predicted by the Paczynski Eq. ([5|) 
and Eggleton Eq. ^ approximations for several bina- 
ries. The table shows that the numerical results are in 
reasonable agreement with these approximations. 

In Figures [5] and [9] we again focus on the binaries 
J1141-6545 and B1855+09 as representatives of the high 
and low WD mass populations. The key result of our 
simulations is that for both of these classes of binaries, 
equilibrium sequences terminate at Roche Lobe overflow, 



TABLE III: The critical separation for Roche lobe overflow 
in units of WD volume radii Rwo for selected binaries. A^ 
is the numerically calculated critical separation and A]^, 
stand for the analytically calculated Ar using Eqs. (Jsj and 
([6]) respectively. 



Object 


4 N 


Al 


A E 


B2303+46 


2.66 


2.75 


2.66 


J1141-6545 


2.79 


2.87 


2.82 


J0621+1002 


3.34 


3.31 


3.31 


J1713+0747 


4.00 


3.91 


3.94 


B1855+09 


4.23 


4.13 


4.16 


J0437-4715 


4.37 


4.29 


4.31 


J0751+1807 


4.94 


4.84 


4.85 


J1012+5307 


4.96 


4.87 


4.87 


B1516+02B 


5.66 


5.58 


5.55 
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FIG. 7: Orbital angular speed Q plotted against binary separation A. Q increases with decreasing A, and is closely approximated 
by Kepler's third law. 




!2M T QM T 

(a)J 1141-6545 (b)B1855+09 

FIG. 8: Binding energy, Ef, according to Eq. (|37|) . rescaled by E^, as defined by Eq. (|36[>. The binary becomes increasingly 
bound with greater angular speed il (smaller binary separations A, see Fig. 0. At very large separations, the energy of the 
system approaches Eoo, as expected. At small separations we see no turning point in equilibrium energy in either case, meaning 
that no ISCO is encountered before reaching the Roche limit. 



rather than at an ISCO or at contact. At the minimum 
separation plotted in the figures our sequences cease to 
converge, and as shown in the contour plots of the fol- 
lowing section, the WD completely fills its Roche lobe. 
We may therefore identify this binary separation with the 
critical binary separation Ar. The graphs of the bind- 
ing energy and the angular momentum do not display a 
turning point, meaning that Roche lobe overflow occurs 
before the binary orbit becomes unstable at an ISCO. 
Also, the minimum separation is larger than the radius 
of the WD and hence the sequences do not terminate at 
contact. 



2. Density contours 



In Fig. [Til] we show the increasing distortion of the WD 
as the binary approaches the critical separation Ar. At 
the critical separation, the WD forms a cusp at the inner 
Lagrange point. Once the binary reaches this critical sep- 
aration, mass transfer from the WD onto the NS across 
this inner Lagrange point will occur. 
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FIG. 9: The angular momentum J versus the angular speed S7. For both of these binaries, J monotonically decreases with 
increasing fl up to the termination of these plots at the Roche limit. This shows that an ISCO is not encountered before 
reaching the Roche limit Ab.- 



1.6 




FIG. 10: Contour plots of the WD in the binary B1855+09 as it approaches the Roche limit. Contours of constant density are 
plotted in the orbital plane. Coordinates are given in grid coordinates, and show the rescaling between iterations such that 
the WD always extends from x = 1 to x = —1 along the rc-axis. In the left panel, at A = 1.23Ar the WD is elliptical, and 
contained well within its Roche lobe. In the middle panel, A — 1.03Ar, the WD is greatly distorted, and nearly fills its Roche 
lobe. We see a sharp cusp forming facing the NS. In the right hand panel, at A — Ar, the WD exactly fills its Roche lobe, and 
has formed a cusp at the inner Lagrange point between the two stars. When the binary separation further decreases the WD 
will overflow its Roche lobe and mass will flow across this Lagrange point onto the NS. 



III. TIDAL DISRUPTION VERSUS STABLE 
MASS-TRANSFER 

As discussed in Section [TlJ GW-driven inspiral drives 
the binary to the Roche limit, at which point mass trans- 
fer from the WD to the NS will commence. There are at 
least two possible scenarios for the subsequent evolution: 

• Mass flows on a secular timescale ~ tow (cf- 
Eq. [TUJ) from the WD toward the NS through the 
inner Lagrangian point while the binary separation 



increases. This process is called stable mass trans- 
fer. 

• The WD gets tidally disrupted by the neutron star, 
leading the system to a merger. This process com- 
prises unstable mass transfer and its timescale is 
of the same order as the orbital timescale of the 
binary, (cf. Eq. [9j. 

Which one of these two scenarios is likely to be realized 
in any particular binary is the subject of the following 
section. 
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A. Stability of mass transfer 

In this section we study the stability of mass transfer 
for a WDNS system near the Roche limit. Our analysis 
follows the approach employed by previous authors, e.g. 
Hut and Paczyhski (54| who studied low-mass semide- 
tached binaries, Verbunt and Rappaport [82], who stud- 
ied WDNS binaries, low-mass X-ray binaries and mass 
transfer from a (sub) giant star, and Faber et al. [53[, 
who studied BHNS binaries. Here we modify the analy- 
sis of Verbunt and Rappaport for WDNS binaries using 
an approximate general relativistic mass-radius relation 
for a cold degenerate WD, rather than the Newtonian 
one which was employed in [82j . 

When the WD component of a WDNS binary loses 
mass, its radius and Roche lobe radius increase. The 
criterion for mass transfer stability is the following: after 
the WD has filled its Roche lobe, instability occurs when 
the timescale for expansion of the WD Roche lobe radius 
is shorter than the timescale for expansion of the WD 
radius, 



flwD Rr 
i?WD Rr 



(44) 



Before we proceed let us distinguish two different mass- 
exchange scenarios [54|, |82j . These are the "no-disk" and 
"with-disk" cases. 

Lubow and Shu (83| demonstrated that once mass 
transfer commences, mass flowing from the secondary 
does not reach the primary directly, but forms a rotating 
Keplerian disk around the primary. According to [54| the 
"no-disk" scenario corresponds to the case where the disk 
is highly viscous and all matter lost from the WD very 
quickly accretes onto the NS. Assuming that the NS is 
not spun up, all angular momentum is then transferred 
very efficiently from the disk back to the orbit. In this 
case a massive disk onto the primary cannot form, hence 
the name of the scenario. This is entirely equivalent to 
the conservative mass-transfer case of Clark and Eard- 
ley [56|. The "with-disk" scenario corresponds to the 
opposite case where the disk is inviscid. As a result an- 
gular momentum cannot be transfered from the disk and 
appreciable mass and angular momentum can be stored 
there forming a massive disk. 

First we calculate each side of Eq. (|4"4")l in the "with- 
disk" case. Assume that mass transfer begins from the 
WD toward the NS at a rate M WD . The WD radius 
depends only on the mass of the WD and hence 



i?WD d In i?WD M> 



WD 



Rwd d In Mwd M 



WD 



(45) 



To account for the effect of the disk on the orbital 
motion we will make the following approximation: We 
assume that the combined gravitational potential and or- 
bital angular momentum of the NS+disk system is the 
same as that of a star of mass M = Mns + Md, where 
Ma is the mass of the disk. In this approximation it is 



the WD with the NS + disk system that is in a Keplerian 
orbit and not just the NS with the WD. Given that the 
disk is gravitationally bound to the NS, this approxima- 
tion is reasonable because the size of the disk is quite 
smaller than the binary separation. In particular, calcu- 
lations by Lubow and Shu (83j show that for mass ratios 
q < 1 the size of the disk edge is only a small fraction 
(2 - 4%) of A. 

We introduce the mass ratio q = A/wd/M. The time 
derivative of q is given by 



— = ( 1 + ^)t? — • 

q AJwd 



(46) 



Using Eq. Q the RHS of Eq. flU]) can be written as 



Rr 



A d\nq q 



(47) 



where we used q instead of q because the WD Roche lobe 
now responds to the combined gravitational potential of 
the NS+disk system. 

We calculate A/ A by imposing angular momentum 
conservation. The disk formed by mass lost from the 
WD, remains in a (equivalent) circular orbit whose ra- 
dius i?d depends on the binary mass ratio q [83[. Under 
these assumptions the total angular momentum of the 
system, J, is given by 



J Jorb ~\~ J disk 3 

where the orbital angular momentum is 

J orb 



M^A^- 

T (l + 9o) 2 



and the disk spin angular momentum is 
Jdisk = (M NS i? d ) 1 / 2 M d . 



(48) 



(49) 



(50) 



Here Mt = Mwd + M and i?d = Ary L1 where is a 
function of the mass ratio q given by [82| 

r h = 0.0883 -0.04858 log q 

+0.1 1489 log 2 q + 0.020475 log 3 q. (51) 

Following Hut and Paczyhski [54[ we neglect the spin 
angular momentum J s of the WD for the purposes of this 
section. Our estimates show that at the Roche limit J s 
is typically a few percent of either J or b or Jdisk because 
J s scales as 



Js 



;( '*wdV 



3 A 
2 A 



i d In -Rwd 
'd\a.Nh 



1 



WD 



WD 



Mwd 

where (3 is a factor less than unity which accounts for the 
central condensation of the star and other factors of the 
binary mass ratio. Note that this last equation holds true 
because we have assumed the binary to be corotational. 
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Using Eq. 
yields 



48]) . angular momentum conservation, J — 0, transfer becomes unstable if 



- = -2 



q-({l+q)r h ) 1/2 



Ah 



WD 



Af\vD ' 



(53) 



where we have dropped terms of order M^/M-^s, 
Ma/MwD (and hence set q Q = q) because at the onset of 
mass transfer M d < M WD and M d < M NS . Eq. Q is 
the same as that of Verbunt and Rappaport [83] (note: 
their q is the inverse of our q). 

Combining Eqs. {H} - (gZJ) and we find that mass 



d In i?wD 
dlnM 



> 2 



WD 



[(1 + <zH] 1/2 



dhxf 

ding 



(1 + ?). 
(54) 

Note that although our derivation so far corresponds 
to the "with-disk" scenario, Eq. l[5"4"]) reduces to the "no- 
disk" scenario, if we set r h = 0. 

To complete the derivation we need to find the loga- 
rithmic derivative of the WD radius , d In Rwu /dhx A/wd ■ 
In [H, HI, HI] the mass-radius relationship for a WD was 
approximated via the Eggleton formula 



J 



i?WD 



Ah 



WD 



AI, 



ch 



-2/3 



M\VD 



2/3-1 1/2 



Ah 



WD 



M„ 



-2/3 



M 



WD 



l-i -2/3 



(55) 



where c = 0.0114, d = 3.5, M c . 



IMMq is the Chan- 



drasekhar mass and M p = O.OOO57Af . In [82| it is stated 
that this approximation corresponds to a fully degener- 
ate WD made of pure He. We have fit Eq. ([5"5]) to the 
mass-radius data plotted in Fig.[5j We find that the WD 
mass-radius relationship is very accurately described by 
Eq. ([55]) for (see Appendix EJ 



c = 0.0116, d = 4.6, 
M c/l = 1.456M , M p = 4.8 • 10~ 5 A/ Q . 



(56) 



By virtue of Eqs. ([53]) and ([55]) we can now determine 
whether mass transfer will be stable or unstable. Follow- 
ing [S2], in Fig. [11] we plot — dlniiwD/dln il%D versus 
the mass ratio of the WDNS binary for several NS masses. 
In addition, we show the RHS of inequality (|54[) calcu- 
lated using the Eggleton approximation to the Roche 
lobe radius. Furthermore, we show both the "no-disk" 
(r/j = 0) and "with-disk" (77, given by Eq. ([51]) ) models. 
The intersection of the — d In i?wD /d In A/wd curve with 
the curves corresponding to the RHS of Eq. ([51]) yields 
the critical mass ratio g cr it- For q > q cr n we have UMT 
(tidal disruption). 

Due to the uncertainty over which of the "no-disk" 
and "with-disk" scenarios actually describes a given bi- 
nary we conclude that there is a minimum and a max- 
imum critical mass ratio, (?crit,min and f? C rit,max, respec- 
tively. When the mass ratio is such that ((?crit,mm < q < 
<7crit,max) a more detailed study is required to determine 
the fate of a binary. Such a study would require precise 
knowledge of the disk viscosity mechanisms. This is a 
task beyond the scope of the current work, but a brief, 
qualitative discussion of this subject can be found in [8^ ]. 
For our purpose it suffices to say that if the mass ratio is 
larger than gcrit.max, then UMT will set in and the WD 
will be tidally disrupted. If the mass ratio is smaller than 
?crit rain) then SMT will set in. 



For a given WDNS binary we can determine (feit.min 
and 9crit,max by numerically solving the following equa- 
tion 



d In -Rwd 
dm M A 



WD 



[{l + q)r h 



1/2 



din/ 
din q 



(! + ?)■ 



(57) 

We solved Eq. ([57]) using a Newton-Raphson algo- 
rithm for 10000 binary models with NS masses in the 
range [0.85, 2.15]M Q . Using these solutions we obtained 
accurate fits for g C rit,min and q C rit,max as functions of the 
NS mass. We find that the following formulae provide 
such excellent fits. 



9crit,max= 0.658 - 0.107M NS 
- 0.034M£ S 



•O.OllM&g, 



gcrit.min = 0.236 - 0.047M N s - 0.003M, 



NS- 



(58) 



(59) 



In Fig. 11(b) we plot formulae ([5"8]) and ([5"9"]) . and the 
numerical solutions of Eq. ([FT]) . We find that the error 
with these formulae is of order 1 part in 10 3 . 

With Eqs. ([58]) and ([59)) at our disposal we can now 
predict the possible fates of observed binaries. We have 
indicated this in our tablesUandHU When the mass ratio 
of a given binary falls in the range feerit.mm, ?crit,max] we 
cannot predict the fate of the binary. 

Finally, an interesting question to ask is what is the 
number of LISA-detectable galactic WDNS binaries that 
will undergo either SMT or tidal disruption per year. 
To answer the question we use the population synthesis 
results of Nelemans et al. [13] as follows: Figure 5 in 
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FIG. 11: (a) Graphs of dlnR/dlnM versus M relation for cold degenerate general relativistic WDs. The open triangles 
correspond to the RHS of Eq. (|54|l in the "no-disk" case. The open circles correspond to the RHS of Eq. (|54p in the 
"with-disk" case. We have used the more accurate Eggleton approximation to the Roche lobe radius for these plots, (b) The 
lines represent qtrit.max and g C rit,min as given by Eqs. (f58|) and ([59)l . The open triangles and circles represent qtrit.max and 
ffcrit.mm obtained by the numerical solution of Eq. (|57p . The insets show the fractional errors of the fitting formulae, where 

Aq = | (/numerical — 9fit | • 



[40J shows the distribution of resolved WDNS binaries vs 
frequency and chirp mass M., where 



M 



(M NS M WD ) 3/5 
(M W d +M N s) 1/E 



= Mi 



7 3/5 



NS 



(1 + g) 1 / 5 ' 



(60) 



We can calculate the binary mass ratio as a function 
of A4, for a given NS mass, via Eq. (|60| . However, Nele- 
mans et al. [40| do not provide the NS masses of the in- 
dividual LISA-detectable WDNS binaries they obtained. 
What they provide is the range of the NS masses which 
is [1.25, 1.55] Mq. They also provide the corresponding 
range of chirp masses, which is [0.4, 1.2]M©. This infor- 
mation is sufficient to place constraints on the number of 
stable or unstable LISA-resolved WDNS binaries. 

We show the results of this calculation in Fig. [12J On 
the q vs M curves we plot the maximum and minimum 
critical mass ratios calculated via Eqs. ([55]) and (|55|) . 
Figure QJ tells us that of the 38 LISA-resolved WDNS 
binaries per year, calculated in [40l| , all those which have 
chirp mass greater than 0.89M Q will undergo tidal dis- 
ruption, whereas those which have chirp mass less than 
0.43M© will evolve into the SMT regime. The fate of 
those in between is uncertain. After assessing the data 
given in we find that after a year of integration LISA 
will resolve about 16 WDNS binaries that will eventually 
undergo tidal disruption and 6 that will undergo SMT. 
However, if the "no-disk" mass-transfer stability crite- 
rion is used, i.e, if we use q C rit,max as the critical value 
for (in)stability, we find that LISA will resolve about 20 
WDNS binaries per year that will undergo SMT. 
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FIG. 12: WDNS critical chirp masses for SMT vs tidal dis- 
ruption. The curves show the mass-ratio as a function of M 
for NS masses that cover the range of NS masses given in 
[40l |. The filled circles and triangles represent the maximum 
and minimum critical mass ratios for SMT, respectively. 



B. Stable mass transfer 

In this section we adapt the approach of Clark and 
Eardley who studied NSNS binary systems, to follow 
the evolution of WDNS systems in the SMT regime. A 
similar approach was also used and adapted for BHNS 
systems in [53|, for double WD systems in [55| and for 
low-mass compact binaries in [57j . 

The Clark and Eardley approach derives an evolution 
equation for the binary by requiring that the orbital an- 
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gular momentum lost from the system be carried away 
by the emission of gravitational waves. However, we will 
also distinguish the case where angular momentum can 
be lost from the orbit and stored in a Lubow-Shu disk 
(the "with-disk" case), in addition to the angular momen- 
tum lost in the form of gravitational radiation. Our basic 
assumptions are: (1) mass is conserved, (2) mass transfer 



takes place across L\ when the WD fills its Roche-lobe, 
and (3) the evolution takes place in quasi-equilibrium 
during which the gravitational radiation reaction can be 
modeled in the quadrupole approximation. 

The sum of orbital angular momentum and the angular 
momentum stored in the disk is given by Eq. (|48|l . The 
time derivative of J is given by 



J 



J 



Jdii 



A_ 

2A 



(1 - do) + ~ 



1 d\nr h Jdisk 



2 dhiq J or b 



— (l + q )rh[q] 



1/2 



Af\vD 



(61) 



where we have used Eq. (|46|) . Note that in contrast to 
[55l ] we retain terms of order Md/AfwD or Md/AfwD- 
These terms are small only near the onset of mass trans- 
fer. Once appreciable mass is stored in the disk these 
terms should not be neglected. 

To complete the calculation we need to link the time 
derivative of the orbital separation to the mass-loss rate 
from the WD. To do this we follow [H[ and solve Eq. |47|) 
for A/ A using the fact that i?wD = Rr at the onset of 
SMT, yielding 



A 
A 



d In i?wD d In / 
d In Mwd d In q Q 



(1 + g ) 



WD 



WD 



(62) 



The angular momentum lost in the form of gravita- 
tional radiation Jew is given by [4l| 



Jgw 32 A/wdMns-Mt 

Jnrb 5 A 



(63) 



Since we assume that mass transfer takes place at the 
point where Roche lobe overflow occurs we have A — 
Ar = Rwr>/f{q ), and hence 



Jgw 



-'orb 



32 (1 + q )f{q ) A Af WD 
R a 2 R 4 



(64) 



Angular momentum conservation implies that 



J 

J orb 



'GW 



(65) 



Substituting Eqs. j6T]), (62|) and (64|) into Eq. ([65]) yields 
the mass transfer evolution equation in the general form 



M W d = - — 



Jd 



d In ]?wd d In / 



d In AfyvD 
+ 2 11 



d In q Q 



(1 + go) 



9o + ■x 



1 d In r h Jd^k 

2 dhiq J or b 

I 



q 



q )r h [q} 



_1 (i + <Zo)/(g ) 4 M 



WD 



qi 



R 4 

"WD 



(66) 



We now need to prescribe the WD mass-radius relation 
and the Roche lobe radius relation. For the former we 
will use Eq. (|55|) and for the latter we will consider the 
more accurate Eggleton approximation Eq. ©, so that 

din/ q^ 3 - q 2 ' 3 + q - 2(1 + g) ln(l + q 1 / 3 ) 



d\nq 



3[c 2 g 2 / 3 +ln(l + «V3)] 



We make the following points: 



(67) 



• Although we derived Eq. l|66|) for the "with-disk" 
case, we can use it to track the evolution of the 
system in the "no-disk" case by simply setting = 



and A/d = = Jdisk- 

• For SMT the initial conditions need to satisfy cri- 
terion (|54p . Equation (|5"^|) arises when the term in 
square brackets in Eq. (|6"6")) is less than zero, for 
Jdisk/ Jorb ^ 1 and q = <?; these conditions hold 
at the onset of mass transfer. Therefore, the closer 
the initial mass ratio is to the critical mass ratio 
for a given binary the faster this binary system will 
evolve. 

• Once the solution for the evolution of the WD mass 
is found we can find the evolution of the separation 
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FIG. 13: Evolution in the SMT regime assuming the "no-disk" scenario for objects J0621+1002, J1713+0747, and B1516+02B. 
(a) Mass evolution of the WD component of the binary, (b) Evolution of the separation of the binary. The curve which 
corresponds to object B1516+02B essentially coincides with the A/A-r = 1 axis during the time span plotted. 



of the binary from Eq. (|6"2")) . When SMT begins 
the separation of the binary will increase because 
the term in square brackets of Eq. (j6"2"|) is always 
negative and hence the RHS of Eq. (f6"2"|) is always 
positive. 

When the separation is known we can also obtain 
the gravitational waveforms by neglecting small de- 
viations from circular Keplerian motion and using 
the quadrupole approximation. For an observer sit- 
uated at r,9,4> = 0, the "plus", h + , and "cross", 
h x , polarizations of the gravitational waves are 
given by 85] 



rh + =-2tt 2 tiA 2 (l+cos 2 9) cos2fit, 



4n 2 fiA 2 cos 9 sin 2flt, 



(68) 



where /i = Mwd^ns/Mt is the reduced mass and 
fl = ^Mt/A 3 the Keplerian angular velocity. 



TABLE IV: WDNS binaries from Tables [TJ [TT] that will un- 
dergo SMT and their expected GW frequencies (in units of 
10 -2 Hz) when Roche lobe overfill takes place and SMT sets 



PSR 


few (Kn 2 Hz) 


B1516+02B 


0.57 


B1855+09 


1.21 


J0437-4715 


1.06 


J1012+5307 


0.70 


J0751+1807 


0.55 


J1232-6501 


0.77 



In Table lTVl we show the expected GW frequencies from 
binaries of Tables Q] and |TT] that will undergo SMT. All 
frequencies are quoted at separation A = Ar. These 
frequencies are above the double WD background noise 
and fall within the LISA frequency range, which is 
[Kr 4 Hz,lHz]. 

We now adopt object B1516+02B as representative 
to study whether these binaries emit detectable gravita- 
tional waves. We obtain a lower limit on the amplitude 
of the waves coming from this object by assuming it is 
located at maximum distance from Earth and 9 = n/2. 
Since these objects are all within the Galaxy we assume 
that the object lies at the opposite end of the Milky way 
with respect to our solar system. The radius of the disk 
of the Galaxy is about 15 kpc and our solar system is 
located 8kpc away from the center of the Galaxy. There- 
fore, the distance from Earth to B1516+02B has to be 
smaller than about 23kpc. If we now use the data of Ta- 
ble IIVI we estimate that the minimum amplitude of the 
expected GW at the Roche limit is /i+, m i n ~ 2 • 10~ 23 . 
LISA's predicted strain sensitivity at 10 -2 Hz is of order 
h ~ 10~ 23 . Thus, WDNS binaries near their Roche limit 
could be detectable by LISA. 

To determine the evolution of a given binary system we 
need to solve Eqs. (|6"2"]l and (|6"6"|) numerically. We adopted 
an adaptive step fourth-order Runge Kutta method to 
carry out the numerical integration. In Fig. [T3]we show 
the results of these calculations, assuming the "no-disk" 
scenario and setting initial conditions corresponding to 
objects J0621+1002, J1713+0747, and B1516+02B from 
Tables UJ and [Til We choose these objects because they 
cover almost the entire range of observed mass ratios for 
SMT in the "no-disk" case. 



In Fig. 13(a) we show the evolution of the WD mass for 
these three systems. We can see that typical timescales 
for the WD to lose half its initial mass range between a 
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few thousands of years (for mass ratios near the critical 
mass ratio) to a few millions of years (for low mass ra- 
tios). As a consequence, typical values of Mwd, i-e. the 
WD mass-loss rate, range roughly between 10 -4 — 10~ 8 
M /yr. If any fraction of this mass accretes onto the NS, 
then it would generate an X-ray photon luminosity that 
would accompany the GW signal from the source. 

In Fig. 13(b) we show both the SMT epoch and the 
earlier inspiral epoch, modeled in the quadrupole approx- 
imation. The evolution equation of the binary separation 
during inspiral is obtained from Eqs. (fBTj) . (|64p and 
setting M\vd = and is given by 



A 3 A = ~y/iM|, 





(69) 



If we integrate Eq. (|69|) and impose the condition that 
A = Ar at t = 0, we find 



A = An 1 



t 



taw 



1/4 



(70) 



Using Eq. (|70[) we can patch the inspiral epoch to the 
SMT epoch. 

Fig. [T3l demonstrates that object J0621+1002 evolves 
on a timescale shorter than that of object J1713+0747 
which in turn evolves on a timescale shorter than that of 
object B1516+02B. We can understand this by checking 
the SMT timescale, which is defined as £smt = A] A, 
given by Eq ([62]) . 

Fig. |13(b)| shows that in the case of object J1713+0747 
when the binary separation is small the inspiral 
timescale, which is defined in Eq. (fT0|) . is longer than 
the SMT timescale, while the reverse is true at large sep- 
arations. To explain this behavior, we study the ratio of 
the inspiral timescale, taw to that of SMT, tsMTi at the 
same separation. 

In Fig. [14] we plot tcw/tsMT as a function of q. 
For high q, t SM T -C t GW , whereas for low q, t SM T > 
tew- The initial conditions which correspond to object 
J1713+0747 yield a q such that at separations close to 
Ar, tew — 2igMT- However, at large separations, after 
the WD has lost substantial mass and q has decreased, 
£gw is smaller than isMT, which explains the behavior of 
the orbital evolution of J1713+0747. Similar reasoning 
explains the behavior of the orbital evolution of object 
B1516+02B. 

Once the numerical solution in the SMT regime is ob- 
tained we can calculate the gravitational waveforms from 
Eq. (j6"5j) . In Fig. \T5\ we show the "plus" and "cross" 
polarizations of the gravitational waves for J0621+1002, 
setting 8 — ir/3. From this figure it is clear that as a 
WDNS binary evolves in the SMT regime both the am- 
plitude and the frequency of the GWs emitted decrease. 
The former results (a) because the reduced mass \x goes 
down, since the mass of the WD decreases with time and 
q < 1, and (b) because the separation increases with 
time. The frequency decreases because of the increasing 
separation. 
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FIG. 14: Inspiral timescale to SMT timescale ratio. 



As mentioned, the numerical results presented so far 
concern the "no-disk", — 0, case. However, the the- 
oretical framework we presented can be applied to the 
"with-disk" case as well. The results when we set ^ 
are qualitatively similar to those of the = case and 
for this reason we will not present them. 

We note that there is one scenario where the quasista- 
tionary treatment of the current section is not applicable 
and instead a fully relativistic simulation is necessary in 
order to follow the evolution of a WDNS system. When 
SMT sets in, mass accretion onto the NS may drive the 
NS above the TOV mass limit. In such a scenario the 
tools of numerical relativity are necessary to follow the 
collapse, which may be delayed, and calculate the emit- 
ted GWs. A very interesting observational aspect of such 
systems is that they offer the exciting scenario that be- 
gins with detectable LISA signals, and ends with LIGO 
signals (from the NS collapse). Note also that object 
B1516+02B is potentially one such system because its 
NS has a mass of 2.O8M . 



C. Unstable mass transfer: Tidal disruption 



As discussed in Section IIII A[ if the mass ratio of a 
given WDNS binary is larger than the critical mass-ratio 
for mass-transfer stability, mass transfer becomes unsta- 
ble. It is important to ask then: what are the possible 
outcomes for this unstable tidal disruption? This is the 
subject of the current section. 

One possible scenario is that the NS plunges into the 
WD and spirals toward the center of the star liberating 
its gravitational potential energy as heat in the WD ma- 
terial. Alternatively the NS may be the receptable of 
massive debris from the disruption of the WD, some of 
which may be shock-heated to characteristic virial tem- 
peratures (kT ~ Mm p /R ~ 10 9 °K). In both cases a 
quasistationary massive NS may be formed with a hot 
envelope. We might assign an adiabatic index of 5/3 to 
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FIG. 15: Gravitational waveforms corresponding to the SMT evolution of object J0621+1002. The curves show the envelope 
of the wave and the insets show the waveforms at different retarded times, (a) The "plus" polarization of the GW, (b) the 
"cross" polarization of the GW. 



the heated gas, i.e., n = 3/2. In this case the EOS no 
longer corresponds to that of strictly cold, extremely de- 
generate electrons (for which n is somewhere between 3/2 
and 3, depending on the mass), but instead corresponds 
to semi-degenerate electrons and thermal ions. We will 
call this the "hot remnant" scenario. 

Notice the resemblance of such an object to another ex- 
otic class of objects, the Thorne-Zytkow Objects (TZOs) 
[86| . These are the end result of a merger between a NS 
and a red giant and might occur in the dense cores of 
globular clusters. The evolutionary path of a TZO in- 
volves mass accretion onto the NS which eventually ex- 
ceeds the OV limit and the core collapses to form a black 
hole. 

Once the hot remnant cools, its fate will depend on its 
mass and spin as well as the angular momentum profile. 
The remnant will settle down into a stable equilibrium 
configuration if its total rest mass is less than that of 
a stable NS with the same angular momentum. Rem- 
nants with rest masses exceeding that of "supramassivc" 
NSs [70j ultimately collapse, unless they are supported 
by differential rotation ( "hypermassive" NSs, [7l|, l72j|). 
Analyzing its fate requires a simulation in full GR, since 
(1) it is only an equilibrium configuration constructed 
in relativistic gravitation that exhibits a maximum mass 
and (2) tracking the evolution requires a dynamical sim- 
ulation. We intend to perform such a simulation and this 
will be the subject of a future work. 

As a first step towards determining the properties of 
these remnants, we will model them as equilibrium con- 
figurations consisting of a point mass NS sitting at the 
center of an extended envelope formed from the WD 
debris. In the "cold" remnant scenario we shall as- 
sume mass and angular momentum (but not energy) 
conservation and employ the ellipsoidal approximation 
(cf. [HI, H3) Hl| ) to treat rotation. Furthermore, we will 



assume that turbulent viscosity (or magnetic fields) act 
to drive the configuration to uniform rotation. For the 
"hot" remnant, we assume that cooling is slow so that 
energy is conserved. 



1. Equilibrium configurations of a point mass surrounded 
by a rotating polytropic envelope 

To construct approximate Newtonian equilibrium con- 
figurations of rotating WDNS remnants which have a 
point mass NS at their center and an extended rotating 
envelope made of the WD debris, we closely follow the 
approach of [41, 88] in which one defines an energy func- 
tional (the total energy of a configuration) and then es- 
tablishes equilibria by extremizing this energy functional. 
The contributions to the energy are (a) the internal en- 
ergy of the envelope, (b) the spin kinetic energy of the 
envelope, (c) the gravitational self-energy of the envelope, 
and (d) the gravitational interaction energy of the enve- 
lope with the point mass at its center. We now calculate 
all these components separately. 

We make three main assumptions: (a) In our energy 
functional we adopt a polytropic density profile param- 
eterized by its central value, as an approximate "trial 
function" for the matter in the envelope; (b) the enve- 
lope is formed entirely from material of the progenitor 
WD and we set its mass is equal to Mwd; and (c) the 
remnant can be modeled as an oblate spheroid having as 
principal axes of its outer surface a\ {— 02) and 03, with 
(J3 measured along the rotation axis and a\ in the equa- 
torial plane. The polytropic density profile is assumed to 
be constant on self-similar spheroids. 

Minimizing the energy functional with respect to cen- 
tral density and eccentricity for fixed mass and angular 
momentum uniquely determines the oblateness and den- 



19 



sity profile of the spheroidal envelope in dynamical equi- 
librium Hi Hi]. 

The total internal energy U of the polytropic envelope 

is Hi 



U = k x Kp x J n M^ 



WD, 



(71) 



where K is the polytropic gas constant and n the poly- 
tropic index, M\vd the mass of the envelope, p c its central 
density and 



(72) 



where £ and 6 are the usual Lane-Emden variables for a 
polytropc and £i and 6± their values at the stellar surface. 
The potential self-energy W of the envelope is [881 ] 



where 



and where 



W 



fc 2 = 



-k 2 M^pl/ 3 g(e), 
3 /4tt|0^ 1/3 



sin 



-(i-e 2 ) 1/a 



9(e) 

is a function of the eccentricity e, denned as 



(73) 



(74) 



(75) 



(76) 



The potential self-energy of the envelope can also be 
written as 188] 



W = -- 



4d 



5-n i? 



9(A). 



(77) 



Here i? is the volume radius of the spheroidal envelope 
given by 

R = ( ai a 3 ) 1/3 , (78) 
so that the central density becomes 

£iM WD 

We have also introduced the oblateness parameter A 

2/3 



A 



= (1-0 



2sl/3 



(80) 



in terms of which the function g becomes 

,g(A) = A 1 /2(i _ A^-^cos-^A 3 / 2 ). (81) 
The rotational kinetic energy T of the envelope is [88] 
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k 3 XJ 2 M, 



-5/3 2/3 
WD rc i 



(82) 



where 



and 



5 (4tt) 



2/3 



4 K« 



2/3 



_ 5 f_l e-e^ 

The moment of inertia I is given by 



I = -K n XM W uR ■ 
5 



(83) 



(84) 



(85) 



If <fi is the gravitational potential due to the envelope, 
the gravitational interaction energy W% between the en- 
velope and the point mass Mns, placed at the center of 
the envelope, is 

Wi = J p(x)0(x)d 3 x = J Af NS< 5(x)0(x)d 3 a; 

NS NS 

(86) 

= 0(O)M NS , 

where 4>(0) is the central gravitational potential due to 

For the ellipsoidal approximation one computes first 
the gravitational potential energy for a spherical (nonro- 
tating) configuration and then "corrects" this energy for 
rotation by multiplying by a correction function. The 
correction function is obtained by considering the rotat- 
ing and nonrotating incompressible case (constant den- 
sity everywhere). We calculate the correction function 
first. 

The gravitational interaction energy W^' 1 between an 
incompressible (constant density p), nonrotating enve- 
lope of mass Mvvd arid radius R, and a point mass M^s 
at its center is 

Wr = -M NS Jf = -l^f™. (87) 

The interaction energy W[' 1 between an incompressible 
rotating envelope of mass M-yvD and volume radius R 
and a point mass Mns at its center is given by Eq. (86) . 
where the potential r,I (O) at the center of the envelope 
is given by [411 ] 



c^(0) = -2tt paf 



sin e, 



where p = 3M WD /(47ri? 3 ). Combining Eqs. (76), (78]), 
(PI). ([SB] and (55) we find 



W r - 1 - 3 M wpM NS sin 1 e _ 2 1/6 
2 i? e 1 j 



3 M WD M] 



(89) 



NS 



1? 



9(A). 



Comparison of Eqs. (57) and (55) yields the rotation 
correction function which is 5(A). 
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In order to find the interaction energy, Wi, between 
a rotating polytropic envelope and a point mass at its 
center, all we are left to do is find the interaction en- 
ergy, W? , between a spherical polytropic envelope and 
the same point mass located at its center, and then ac- 
count for rotation by multiplying by g{\). Given Eq. 
(155)1 . this implies that we only have to calculate the po- 
tential (f)(0) of a polytrope at its center. 

The gravitational potential at a point x of any arbi- 
trary matter distribution p(x) is given by 



Therefore, at x = the potential is given by 

dm 



0(0) = - 



(90) 



(91) 



where dm — 4irpr 2 dr. Integration of Eq. (|91|) by parts 
yields 



0(0) 



Ah 



WD 



-dr. 



(92) 



By using the hydrostatic equilibrium equation we rewrite 
Eq. (H2D as 



0(0) 



Ah 



WD 



R 



1 dP , 
— —dr. 
p dr 



(93) 



Substituting Eq. (|I7|) in Eq. (J93J) we obtain 
Mwd 



0(0) 



R 
Mwd 
R 



K i/T I P -UrdP dr 
dr 



- (n + l)Kp, 



1/r 



where p c is the density at the center and T = 1 
Using the familiar polytropic relations 



R = 

and 

M W d = 47r 
we obtain 



(n + l)K 



4tt 



1/2 



,(l-n)/ 2 n 6i 



(n+ l)K 



-in 



3/2 



WD 



5 (3-„)/2n^2|^| 



(n + l^pVn^i^i 



R 



Substituting Eq. ([97]) into Eq. J94} then yields 



0(0) 



M- 



WD 



R 



1 



1 



(94) 
1/n. 
(95) 

(96) 
(97) 
(98) 



Finally, substituting Eq. (|98|) into Eq. (|86|) we obtain the 
interaction energy, VF/, in the spherical case 



Wf = - 



1 



1 



(99) 



As an immediate check we can insert £i = y/6 and |^| = 
V6/3 for the incompressible case into Eq. (|M)l . which 
yields, as expected, Eq. ([87)1 . 

The expression for the interaction energy between the 
spheroidal polytropic envelope and the point mass at its 
center is obtained by correcting Eq. (|9"9"|) for rotation, so 
that 



IT"; 



1 



M NS M WD 

^4 , r 5/3 1/3 



5(A) 



(100) 



where we used Eq. (f75|) in the last step, and where 



1 + 



1/3 



(101) 



and, as before, q = Mwd/Mns- 

We now write the total energy functional of the point 
mass and spheroidal polytropic envelope as 



(102) 



E{M WB ,g, J;p c ,A) = U + W + T + Wi, 



where U,W,T and W, are given by Eqs. ([7T]), fH?]). 
and (|100p respectively. Given M W d, 9 an< 4 the only 
parameters required to describe an equilibrium configu- 
ration are the central density of the envelope p c and the 
oblateness parameter A. The equilibrium conditions then 
are 



dE n a dE n 
— = and -^r = 0, 

Op c OA 



(103) 



where the partial derivatives are taken keeping q, Mwd 
and J constant. 

The first equality in Eq. (| 103[) yields the virial relation 



n 



2T + W t = 0. 



Here the total gravitational potential energy Wt 

T^;is 



(104) 
: W + 



W t 



where 



<1 



id 

q 



M 2 



R 



1 



5-n q\ ^\9[ 
The second equality in Eq. (|103|) yields 



|Wtl 



1 



3A 3 



3A 3 / 2 



1-A 3 (1- A 3 ) 1 /2 C0S -i A 3 / 2 



(105) 
(106) 

(107) 



which is the usual relation between T/|W| and eccentric- 
ity for a uniformly rotating spheroid (cf. Eq. (7.3.24) in 
[41| ) except that now W is the total gravitational poten- 
tial energy. 
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By virtue of Eq. (|104[) we can obtain the mass-central 
density relationship 



Ah 



WD 



1 



5(A) 1 



2T 



-3/2 



(108) 



where M a is the mass of a nonrotating spherical poly- 
tropic star which has the same polytropic constant K, 
polytropic index n and central density p c as the enve- 
lope. Rewriting Eq. ([96]) in terms of the constants ki 
and lt2, this mass becomes 



\ nk 2 J Pc 



(109) 



Using Eqs. (|79[) and (|108p we can write the equilibrium 
volume radius as 



R — i? 



1 



k2q 



5(A) 1 



2T 



V(3- 



, (110) 



where R is the radius of a nonrotating spherical poly- 
tropic star with the same polytropic constant K, poly- 
tropic index n and mass as the envelope. From Eq 
we can express R Q as 



R Q — b„ 
where 



~(n + l)K~ 


n/(3-n) 




4-7T 




V 47T J 



(l-n)/(3-n) 



-(l-n)/(3-n) 



, (HI) 



(112) 



Finally, if we combine Eqs. (|102|) and (| 1 04[) . the equi- 
librium energy can be written as 



E eq = — -W t \ I 



3 — 2n T 
3-rc Wt 



(H3) 



Eqs. (fT04|) and (fT05)) - pT3)) completely determine the 
equilibrium configuration corresponding to a given Mwd, 
q and J. 



2. Initial energy and angular momentum 

We have assembled all the expressions necessary to 
construct approximate equilibrium configurations of a 
point mass NS surrounded by an extended envelope, and 
will now use these expressions to model the properties of 
the remnants of WDNS mergers. 

We first need to determine the initial energy and an- 
gular momentum of the binary system just before tidal 
disruption sets in at the Roche limit. We treat the WD 
as a corotating, (nearly) spherical polytrope with poly- 
tropic constant K' and index n' (which may be different 
from the K and n of the remnant), and model the NS as 
a point mass. The angular momentum of the binary is 
given by the following formula 



J 



(114) 



where I cm is the moment of inertia of the binary calcu- 
lated in the center of mass frame. We write this moment 
of inertia as the sum of an orbital component I or b and 
the WD's spin component I s , 



with 



and 



'WD' 



1+3 



Is = -K„'MwD-RwD- 

Combining the last two expressions we obtain 
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In Section III El we demonstrated that deviations from 
Keplerian motion are small even at the Roche limit. 
Thus, we approximate the angular velocity of the binary 
by using Kepler's third law 



2 q + 1 M WD 



<l 



A3 ' 



(119) 



Using Eq. Q we can express the critical separation 
A = as 



A, = 



Rwd 

/(?) ' 



(120) 



If we now combine Eqs. (|114p through (|120p we find that 
the initial angular momentum becomes 



j2 _ ^wd%d 



<K? + !)/(<?) 
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The initial energy of the binary is given by Eq. (7.11) 
in [83 as the sum 



E = E p + T' + W[. 



Here E p is the "intrinsic" energy of the polytrope 



EL 



3 - ri 



W 



1 - 




) ^ 1 






I 3 - n> , 


' \W'\_ 



(122) 
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where W is the gravitational self-energy of the WD given 
by 



W 



3 M 2 

JU WD 



5 — n' i?wD 
T s is the spin kinetic energy of the WD 

1 



(124) 



(125) 
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T' is the orbital kinetic energy, 

„ (M wu )A 2 n 2 



T 



2(1 + 9) 



(126) 



and W[ the interaction energy between the point mass 
NS and the polytropic WD 



Wi 



A/wdM N s 
A 



(127) 



where we neglect tidal interaction terms. 
Using Eq. (fTT9|) we can write Eq. (fl~22)) as 



2qA 



(128) 



By virtue of Eq. ()119|) at A = the spin kinetic energy 
becomes 

Ts = ^(i±i)/!M^ (129) 

5 q -Rwd 



Combining Eqs. p24| . (fl~28"f and (fl29|) . we can finally 
write the initial total energy at the Roche limit as 



E = -(l-A) 



3-n'\M WD 



5 — n' J i?wD ' 



(130) 



where 



2(3-2n')(5-n'K,(g + l)/ 3 (g) 
15 ? [2«(3-n') + /( 9 )(5-n')] ' 



3. The "hot" remnant 

To determine the properties of the hot remnant, we 
assume conservation of E, J and Mp = Mwd+Mns. We 
further assume that after tidal disruption all of the WD 
debris forms the envelope surrounding the point mass NS. 
Accordingly, Eqs. ||H2J) and (fTUS")) yield 
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5J 2 A 



\W t \ iK n k 5 M^ D Rg{Xy 



(132) 



where the envelope is assumed to have a polytropic index 
n, constant K and volume radius R. Angular momen- 
tum conservation amounts to substituting J 2 in Eq. (jl32[) 
from Eq. (fT2T]) . We then find 
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WD 

R 



(133) 

We now impose energy conservation. This amounts 
to setting the RHS of Eq. (TTT3j) equal to the RHS 



of Eq. (|128p . Combining the resulting equation with 
Eq. (|105[) we obtain 



i?WD 
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Equations pH7| . p^]) . and (fT34]) form a system of 
three algebraic equations for the three unknowns T/|Wt|, 
A and R/Rwr>- After we specify n' and K' (see below for 
cases considered in this work) we solve this system via 
an iterative Newton-Raphson scheme. Once we obtain 
T/|Wt|, A and i?/i?wD we can determine all properties 
of the hot remnant as follows: Given the mass of the 
envelope Mwd, we find the total potential energy Wt, 
the central density p c and the volume radius of the enve- 
lope from Eqs. (|105|) , (|108j) and (|110p . respectively. Also, 
given T/|Wt| and Wt, we find T and in turn we calculate 
the angular velocity of the rotating remnant as 



(135) 



where / is the remnant's moment of inertia, given by 

Eq. dHHD- 

Furthcrmore, given the initial K' , n' and M^vd we cal- 
culate the final polytropic constant K as follows: The 
radius of the initial WD is given by Eq. (jllip , except for 
constants K' and n' 
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Thus, the ratio of R to i?wD is 
Ro K ( Mwd^W^W / 1 



-Rwd K> \ 47r 
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In the limiting case n = n' Eq. (|137[) reduces to 
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(5(n, n'), 
(137) 

(138) 
(139) 

(140) 



Therefore, once we determine R/Ry/D, Eqs. (|110|) and 
(|137p form one equation which we use to determine K 
(which is related to the final entropy). 
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TABLE V: WDNS binaries from Tables U and [TT] which will undergo UMT. The table shows the properties of the initial WD 
and those of the hot remnant, for which we always assign n = 1.5. The initial WD is described as an n' = 1.5 (low-mass) 
or an n' = 2.9 (high-mass) polytrope. The columns from left to right give the name of the object, the central density of the 
initial WD, the radius of the initial WD, the central density of the envelope of the remnant, the ratio of the final (volume) 
radius to the initial WD radius, the ratio of the spin kinetic energy to the total gravitational potential energy of the remnant, 
the dimensionless angular momentum of the remnant (in geometrized units) , the spin angular velocity of the remnant and the 
central temperature of the envelope. 









n' = n — 


1.5 










PSR 


p?(10 B g/cm 3 ) 


-RwD(Km) 


pi (10 5 g/cm 3 ) 


R/ Rwo 


T/\W t \ 


J/Mt 




9 C (10 9 K) 


B2303+46 


6.88 


8136 


5.11 


2.38 


0.26 


19.7 


9.22 


1.46 


J1157-5114 


5.29 


8500 


3.26 


2.53 


0.27 


20.8 


7.74 


1.26 


J1141-6545 


4.23 


8821 


2.44 


2.59 


0.27 


22.0 


6.82 


1.10 


J1435-60 


4.92 


8602 


2.76 


2.61 


0.27 


20.8 


7.30 


1.23 


J1453-58 


4.66 


8682 


2.51 


2.65 


0.27 


21.0 


7.04 


1.19 


J1022+1001 


3.09 


9295 


1.26 


2.91 


0.28 


22.2 


5.36 


0.96 


B0655+64 


2.70 


9510 


0.99 


3.01 


0.29 


22.5 


4.89 


0.90 








ri = 2.9, n 


= 1.5 










PSR 


p°(W 8 g/cm) 


-Rwd(Khi) 


pt (10 5 g/cm) 


R/Rwd 


T/\Wt\ 


J/M% 


Q( S - 1 ) 


e c (10 9 K) 


B2303+46 


3.38 


4398 


1.61 


6.53 


0.12 


14.0 


4.52 


1.67 


J1157-5114 


0.68 


7202 


0.31 


6.62 


0.13 


18.6 


2.14 


0.97 


J1435-60 


0.54 


7685 


0.24 


6.67 


0.14 


19.1 


1.96 


0.91 



Finally, given K we also obtain an estimate for the 
temperature of the envelope O, if we approximate the 
total pressure P as the sum of the cold pressure P' and 
thermal pressure Pth- This gives the correct form in the 
extreme cold (degenerate) and extreme hot (Maxwell- 
Boltzmann) limits. Then the thermal pressure is 

^ = - - — = K p l ' n - K'p l ' n ' . (141) 
P P P 

The temperature is then given by 



k p 

Here m u is the atomic mass unit and p, is the mean molec- 
ular weight, which, for fully ionized material, has the 
value p, — 12/7 and k Boltzmann's constant. As noted at 
the beginning of Section IHI C[ we will always assume an 
adiabatic index T = 5/3 (n = 3/2) for the hot remnant. 

In Table fVl we list the properties of WDNS remnants in 
the hot scenario by modeling those binaries from Tables 
U and [TT] which will eventually undergo UMT. We have 
considered two different models: (a) The initial WD is 
described by n' = 1.5, and (b) the initial WD is described 
by n' = 2.9. 

For case (a) we set K' = 3.161 x 10 12 (c<?s), which corre- 
sponds to the nonrelativistic limit of the EOS for an ideal 
degenerate electron gas. The data in the table shows that 
the volume radius of the hot remnant is always larger 
than the radius of the progenitor WD, i.e. the hot rem- 
nant is puffed up. The spin frequency of the remnant is 
of order few Hz and J/M^ is of order 20 and hence sig- 
nificantly larger than unity. This result implies that the 



remnant cannot collapse to form a Kerr black hole unless 
some process either removes angular momentum from the 
system, or forms a massive disk about a hole. Also, the 
ratio T/|Wt| is almost always comparable to 0.27. This 
implies that these configurations are marginally prone to 
the development of the bar mode instability on a dynam- 
ical timescale [H [M HI ■ 

The initial and final densities we calculated are consis- 
tent with the choice of polytropic EOS we made, because 
they are not much larger than 10 6 g/cm 3 , below which a 
n = 1.5 index describes well the cold degenerate matter 
[4lJ. However, note that the WD components of objects 
B2303+46, J1157-5114 and J1435-60 are large enough 
(> 1.1M©) that the degenerate electrons may be rela- 
tivistic. If this is so, a polytropic index n' ~ 3 better 
describes the WD EOS. 

To check this we used our p c - Mwd data from the 
integration of the TOV equations in conjunction with 
the EOS, described in Section III B 21 Wc find that the 
central densities which correspond to the masses of the 
WD component of objects B2303+46, J1157-5114 and 
J1435-60 are 3.38 x 10 8 g/cm 3 , 6.75 x 10 7 g/cm 3 and 
5.36 x 10 7 g/cm 3 , respectively. These densities are sig- 
nificantly larger than 10 6 g/cm , and hence these high 
mass WDs arc best studied by a polytropic index n' ~ 3. 

We choose n' = 2.9 to model these quasi-relativistic de- 
generate cases. This value is close to that of an extreme 
relativistic ideal degenerate electron gas and at the same 
time avoids singularities that arise for n' = 3. Also, to 
have a polytropic EOS P = K' p( 1+1 / n ) which is consis- 
tent with p c and Mwd, we choose K' as follows: Usin^ 
the central densities 3.38 x 10 8 g/cm , 6.75 x 10 7 g/cm 
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TABLE VI: The cold scenario for WDNS binaries of Table IVl We show the properties of the cold remnant for ri — n — 1.5 
and ri — n — 2.9 (for the highest WD masses only). The columns from left to right give the name of the object, the central 
density of the envelope of the remnant, the ratio of the final (volume) radius to the initial WD radius, the ratio of the spin 
kinetic energy to the total gravitational potential energy of the remnant and the spin angular velocity of the remnant. The 
central density and radius of the initial WD and the dimensionless angular momentum of the remnant are the same as those 
listed in Table |V] 



ri — n = 1.5 


PSR 


p% (10 6 g/cm 3 ) 


R/ Rwd 


T/\W t \ 




B2303+46 


3.72 


1.23 


0.36 


21.8 


J1157-5114 


2.60 


1.27 


0.38 


18.5 


J1141-6545 


2.01 


1.28 


0.38 


16.4 


J1435-60 


2.31 


1.29 


0.38 


17.6 


J1453-58 


2.14 


1.30 


0.39 


17.0 


J1022+1001 


1.25 


1.35 


0.40 


13.2 


B0655+64 


1.04 


1.37 


0.41 


12.1 


ri = n = 2.9 


PSR 


pi (10 7 g/cm 3 ) 


R/ Rwd 


T/\W t \ 


fi^ 1 ) 


B2303+46 


5.83 


1.80 


0.35 


66.9 


J1157-5114 


1.01 


1.89 


0.36 


28.5 


J1435-60 


0.74 


1.93 


0.37 


24.9 



and 5.36 x 10 7 g/cm 3 of the WDs and the correspond- 
ing masses we compute a consistent K' for each of these 
objects separately, using Eq. and setting ri = 2.9. 

The last three rows of Table fVl correspond to ri = 2.9. 
The volume radius of the hot remnant is again larger than 
the radius of the progenitor WD, but in this case the hot 
remnant is more puffed up than in the ri = n = 1.5 case. 
The spin frequency of the remnant is of order few Hz, but 
less than the corresponding n' = n = 1.5 case. Further- 
more, J/M^ remains much larger than unity, hence such 
remnants cannot collapse to form a Kerr black hole, un- 
less some process removes angular momentum or forms a 
massive disk about a hole. In contrast to the ri = n = 1.5 
case, the ratio T/|Wt| is less than 0.27, but close to the 
secular instability limit for bar formation [§3, [H, [9(| • 



4- The "cold" remnant 

We determine the properties of cold remnants as fol- 
lows: First, we assume that angular momentum is con- 
served and that the hot remnant cools by radiating the 
excess energy until K' = K and n' = n, thereby return- 
ing to its original (degenerate) state. Equations (|107p . 
(jllOp and (|133[) form a system of three equations for the 
three unknowns T/|Wt|, A and i?/i?wD- This is so, be- 
cause in Eq. (|110p R D = i?wD, which follows from Eq. 
(|137[) under the assumption that K' — K and ri — n. We 
solve this system of algebraic equations via an iterative 
Newton- Raphson scheme. Once we obtain T/|W t |, A and 
i?/i?WD 7 we determine all properties of the hot remnant 
as follows: Given the mass of the envelope AfyvD and 
n, using Eqs. (|105|) . (|108p we calculate the total poten- 



tial energy Wt and the central density p c of the envelope 
respectively. Finally, given T/|W t | and W t , we find T 
and in turn calculate the angular velocity of the rotating 
remnant via Eq. (| 1 35[) . 

In Table I VII we list the properties of WDNS remnants 
in the cold scenario by considering the same binary sys- 
tems as those listed in Table [V] Just like in the hot 
scenario we consider two cases (a) n' = n = 1.5 and (b) 
n' =n = 2.9. 

For ri = 1.5 we again set K' = K = 3.161 • 10 12 (c 5 s). 
Note that the central densities are consistent with the 
choice of polytropic EOS we made. The data in the ta- 
ble shows that the volume radius of the cold remnant is 
always slightly larger than the radius of the progenitor 
WD, i.e. the remnant is slightly puffed up. The spin 
frequency of the remnant is of order of a few tens of Hz 
and J/M| is of order 20. Therefore, these objects, too, 
cannot collapse to form a Kerr black hole unless some 
process removes angular momentum from the system or 
forms a massive disk about a hole. The ratio T/\Wt | is al- 
ways significantly larger than 0.27, indicating that these 
cold configurations are prone to the development of bars 
on a dynamical timescale. 

Just like in the hot case, for the three highest mass 
WDs we also use n' = 2.9 and K' = K calculated as 
described in the previous section. We list the results of 
this analysis in the last three rows of Table [VTl Note that 
the central densities of the cold envelope are marginally 
large enough to justify the choice n = 2.9 perhaps with 
the exception of object J1435— 60. The data in the table 
shows that the volume radii of the cold remnants are a 
little larger than the corresponding radii in the ri = n = 
1.5 case. The spin frequency is again of order of a few tens 
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of Hz and larger than the corresponding values of Table 
IVIl These objects have precisely the same J/M| as in 
the n' = n = 1.5 case and hence, the remnants cannot 
collapse to form a Kerr black hole unless some process 
removes angular momentum from the system or forms a 
massive disk about a hole. Finally, the ratio T/|Wt| is 
not very different than the n' — n = 1.5 case. Thus, all 
these cold configurations are prone to the development 
of the dynamical bar mode instability. 

We emphasize again that all these estimates are tenta- 
tive, and meant only as speculative previews of the pos- 
sible characteristics of a WDNS remnant. Clearly, fully 
relativistic, dynamical simulations are needed to explore 
these scenarios conclusively. 



5. Alternative outcomes 

So far we have considered the possibility that the evo- 
lution of an unstable WDNS binary results in an equilib- 
rium remnant with a NS at its center. However, this is 
not the only possible outcome. Alternatives include: 

• Prompt collapse to a black hole immediately fol- 
lowing NS merger with the disrupted WD; 

• Formation of a differentially rotating, hypermassive 
NS following merger that ultimately undergoes de- 
layed collapse to a BH; 

• Ejection of appreciable WD debris, leaving a rotat- 
ing NS remnant. 

We point out that the situation here is reminiscent 
of the situation found in binary neutron star mergers, 
where either prompt collapse to a BH or the formation 
of a hypermassive star followed by delayed collapse occurs 

The frequency of the gravitational radiation during the 
NS collapse phase is of order \JGp (where p is the average 
NS density). If a bar is formed as the NS collapses the 
frequency of the GWs is of the same order as the spin fre- 
quency of the collapsing NS [9l|. Hence, the GW signal 
from a collapsing NS could be between a few tens of Hz 
to a few kHz, suitable for detection by Advanced L1GO. 
Therefore, just like in the UMT case with NS collapse, we 
could again have a scenario that begins with detectable 
LISA signals (from the WDNS inspiral), and ends with 
either a burst or quasiperiodic signal of GWs (from the 
NS collapse) detectable by ground-based interferometers. 
Note that the UMT case would have a different GW sig- 
nature than that of the SMT case. Therefore, these two 
scenarios are not degenerate and gravitational wave ob- 
servations alone could discern which of the two scenarios 
took place. 

Another interesting aspect of unstable WDNS bina- 
ries is that they might also be progenitors of gamma-ray 
bursts (GRBs). The material from the tidally disrupted 
WD may eventually form an accretion disk onto a BH 



remnant, if the NS collapses, whose lifetime could be of 
order typical GRB timescales. 

On the other hand, if the lifetime of the disk is long 
enough and the conditions favorable for fragmentation to 
occur, WDNS mergers may offer a plausible route to the 
formation of planetary systems around isolated neutron 
stars. Currently, one such planetary system is known [92| 
and the formation process of these systems remains an 
open question. 

Fully relativistic hydrodynamics simulations are nec- 
essary to resolve the final fate of WDNS remnants. Al- 
though the field of numerical relativity has now matured 
enough to be able to handle problems involving BHBH 
and BHNS binaries, studying WDNS binary system of- 
fers, yet, an extra complication. This particular problem 
is extremely challenging because it involves a NS of char- 
acteristic size ~ 10 km and dynamical timescale ~ 1 ms 
and a WD of size ~ 10 3 km and dynamical timescale ~ 1 
s. Therefore, in order to be able to resolve all the different 
scales the full power of numerical relativity in conjunc- 
tion with adaptive mesh refinement has to be employed. 
This is what we are now planning. 



IV. SUMMARY AND DISCUSSION 

In this paper we considered WDNS binaries and "set 
the stage" for fully relativistic hydrodynamic simulations 
of the WDNS merger. Like NSNS binaries, but unlike 
BHNS and BHBH binaries, WDNS binaries are known 
to exist and are abundant. We have compiled a list of 
observed WDNS binaries and their properties in Tables 
Hand OH 

The emission of gravitational radiation will cause the 
orbital separation of a WDNS binary to shrink. We mod- 
eled corotational WDNS binaries in circular orbit and 
found that these models terminate at the Roche limit. 
At this point the binary can undergo either SMT and 
evolve on a secular timescale or UMT and evolve on a 
hydrodynamical timescale. 

Following the stability analysis of Verbunt and Rap- 
paport [82j, we showed that the subsequent fate of the 
binary is determined by a critical mass ratio. Using the 
results of this stability analysis we predicted the possible 
fates of known WDNS binaries and we indicated this in 
our Tables Q] and [TTJ Furthermore, based on population 
synthesis results by Nelemans et al. [Ioj|, we gave an es- 
timate of the number of LISA-resolved galactic binaries 
per year that will undergo SMT and the number of those 
that will undergo tidal disruption. We found that ap- 
proximately up to 16 WDNS binaries will undergo UMT 
(tidal disruption), and up to 20 SMT. 

In the case of SMT, the timescale of the mass transfer, 
and hence the timescale of the binary evolution, is set by 
the emission of gravitational radiation. We treated this 
quasistationary SMT epoch by applying the approach of 
Clark and Eardley [56| and Faber et al. [53j , also adopted 
in [HH, [I?], Hll , and we estimated the corresponding grav- 
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itational waveforms. 

We also constructed approximate equilibrium config- 
urations of rigidly rotating WDNS merged remnants in 
order to explore possible outcomes of WDNS mergers. 
We found that unless some process removes angular mo- 
mentum from the system or leads to the formation of 
a massive disk, these massive equilibrium configurations 
cannot collapse directly to form a Kerr black hole since 
their angular momentum ( J/M 2 ~ 20 > 1) exceeds the 
Kerr limit. Furthermore, in most of our case studies we 
found that the merged remnants have a ratio of spin ki- 
netic to gravitational potential energy greater than 0.27, 
which implies that they are dynamically unstable to bar 
formation. However, we emphasize that these prelimi- 
nary results are at best approximate and must be con- 
firmed and/or revised by detailed numerical simulations. 

The fate of a merged WDNS binary depends on the ini- 
tial mass of the cold progenitor stars, the degree of mass 
and angular momentum loss during the WD disruption 
and binary merger phases, the angular momentum pro- 
file of the WDNS remnant and the extent to which the 
remnant gas is heated by shocks as it pours onto the NS 
and forms an extended, massive mantle. These issues all 
require a hydrodynamic simulation to resolve. 

Moreover, ascertaining whether or not the NS ulti- 
mately undergoes a catastrophic collapse to a BH (either 
prompt or delayed) requires that such a simulation be 
performed in full general relativity. In fact, even the fi- 
nal fate of the NS in the alternative scenario in which 
there is a long epoch of SMT may also lead to catas- 
trophic collapse, if the final NS mass exceeds the TOV 
mass limit. This scenario too will require a general rel- 
ativistic hydrodynamic simulation to track. We plan to 
explore some of these alternative scenarios in detail in 
the future, aided by simulations that employ our AMR 



relativistic hydrodynamics code. 



Acknowledgments 

It is a pleasure to thank R. Webbink, V. Kalogera and 
B. Willems for helpful discussions. MM gratefully ac- 
knowledges support from the Maine Space Grant Consor- 
tium. This work was supported in part by NSF Grants 
PHY02-05155 and PHY06-50377 as well as NASA Grants 
NNG04GK54G and NNX07AG96G to the University of 
Illinois at Urbana-Champaign, and NSF Grant PHY07- 
56514 to Bowdoin College. 



APPENDIX A: ACCURACY OF THE 
APPROXIMATE WD MASS-RADIUS RELATION 

In this appendix we demonstrate that Eq. (|55[) (com- 
bined with Eq. (|56p) is a good fit to the actual WD mass- 
radius data we obtained after numerically integrating the 
TOV equations adopting the EOS described in Section 
QTBjfor [i e = 2. 

We find that the fractional error in the radius is less 
than 2% for M W d <E (0.02, 1.39)M© and the fractional 
error in the logarithmic derivative e?lni?/dlnM is less 
than 4.85% for A/ W d € (0.02, 1.33)M S (see insets in 
Fig. Q1)J where we show the fractional errors) . Note that 
the mass range over which the formula (|55l) (combined 
with Eq. (p)6"|) ) is accurate are those which are of astro- 
physical interest for WDNS binaries. We also note our 
modified expression for the relativistic mass-radius rela- 
tion represents an improvement over the original Eggle- 
ton approximation. 
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